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ABSTRACT 

We employ a multipole-Taylor expansion to investigate how tightly the gravitational 
potential of the quadruple-image lens MG J0414-I-0534 is constrained by recent VLBI 
observations. These observations revealed that each of the four images of the 
background radio source contains four distinct components, thereby providing more 
numerous and more precise constraints on the lens potential than were previously 
available. We expand the two-dimensional lens potential using multipoles for the 
angular coordinate and a modified Taylor series for the radial coordinate. After 
discussing the physical significance of each term, we compute models of MG J0414+0534 
using only VLBI positions as constraints. The best-fit model has both interior and 
exterior quadrupole moments as well as exterior m = 3 and m = 4 multipole moments. 
The deflector centroid in the models matches the optical galaxy position, and the 
quadrupoles are aligned with the optical isophotes. The radial distribution of mass 
could not be well constrained. We discuss the implications of these models for 
the deflector mass distribution and for the predicted time delays between lensed 
components. 



Subject headings: gravitational lensing — methods: data analysis — galaxies: structure 



1. Introduction 



There are now about 40 observed examples of "strong" gravitational lensing, in which multiple 
images of a background object are produced by the gravitational potential of an intervening object 
( Kochanek et al. 1999| ). Constructing models for the gravitational potential of the deflecting mass 
in these systems is of fundamental importance, in the first place to verify that the observed system 
is truly a gravitational lens. A crucial hurdle for every promising gravitational lens candidate 
is to have its image configuration reproduced, at least approximately, by a plausible model for 
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the deflecting mass distribution. If this cannot be achieved, the lensing interpretation must be 
seriously doubted. 

Once this and other tests for lensing are passed successfully, however, model construction 
moves beyond a plausibility check into a direct measurement of the mass distribution of the 
lens. Lensing makes a unique contribution because it is sensitive to all forms of mass (including 
dark matter), yet does not depend on any luminous tracers in the lens. Lens modeling is also 
a crucial step in the enterprise of determining cosmological parameters by measuring the time 
delays between the light curves of multiple images. A successful model will predict the values of 
these time delays in terms of parameters such as Hq, 0,rn, and ^}\. These parameters can then be 
constrained by the time-delay measurements ( [Refsdal 1964, 1966 ; for recent examples see Haarsma 



et al. 1999, Lovell et al. 1998|). The appeal of this technique is that it does not make use of the 



usual chain of intermediate distance indicators and their associated uncertainties. 

The determinations of mass distributions and cosmological parameters are both frustrated 
by the main challenge of lens modeling: it is a poorly-constrained inverse problem with unknown 
systematic errors. It is neither obvious which parameterized model for the gravitational potential 
should be chosen, nor how precisely the parameters are constrained by the data. The choice of 
parameterization for the gravitational potential is usually dictated by two factors: preconceptions 
about the mass distribution of galaxies (based on the distribution of luminous matter and/or 
velocity dispersions), and ease of computation. Some examples are singular isothermal spheres, 
elliptical density profiles, elliptical potential contours, and triaxial density profiles (see e.g. 
[Schneider fc Weiss 199T| ; |Nair fc Garrett 19971 ; |Chae, Khersonksy fc Turnshek 1998|) . Added 



to these are terms representing external "shear" from the tidal forces of neighboring mass 
concentrations. Such models have all been used successfully to reproduce the image configurations 
of the known sample of lenses, at least qualitatively, which leads one to wonder how well the 
observations of strong lensing actually constrain the lens potential. 

Kochanek (1991) was the first to investigate this question systematically. He attempted to 
model each of 12 different lenses using a set of 6 different parameterized potentials. The potentials 
he considered were point masses and singular isothermal spheres, added to shear terms with one 
of three radial dependencies. He found that the total mass in the region between the multiple 
images is well constrained, but the radial dependence of the potential is not. This was because the 
multiple images are typically located close to the "Einstein ring" of the lens, so the constraints on 
the potential are correspondingly limited to a small range in radius. Kochanek suggested (but did 
not carry out) a Taylor expansion, parameterized by the deviation from the ring radius, as a way 
to turn this fact to advantage. 

The purpose of this paper is to elaborate upon this suggestion and apply it to a particular 
quadruple-image gravitational lens, MG J0414-I-0534. We expand the potential as a series of 
multipoles in angle, and as a modified Taylor series in radius. The resulting series has three 
advantages over traditional lens models. One, it is mathematically general, and therefore less 
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subject to preconceptions about the galactic mass distribution (which, after all, is what is trying 
to be determined). Two, the degeneracies in the lens equation and the moments of the mass 
distribution have simple correspondences to individual terms of the series. Three, each successive 
term in the series is expected to contribute a smaller amount to the light deflection, so long as 
the radial parameter is small and the angular dependence of the potential is smooth. This allows 
the complexity of the model to be easily prescribed by the truncation of the series. However, 
although the first assumption (the smallness of the radial parameter) is empirically true for 
quadruple-image lenses, the second assumption (the angular smoothness of the potential) is open 
to question. Nevertheless, in the face of our ignorance of true galaxy structure, we employ a 
multipole expansion because of its simplicity and generality. 

The disadvantage of a mathematically general expansion (rather than one that is motivated 
by astrophysical preconceptions) is that the number of parameters is large. In many cases, 
especially for the double-image lenses, the number of parameters would be comparable to or more 
than the number of observational constraints. This is probably why such an expansion has not 
been used previously. However, recent VLBI images of MG J0414-I-0534 provide a much larger 
body of constraints than were previously available ( [Trotter 1998 ). Each of the four lensed images 



was found to contain four components, with clear correspondences between the components of 
different images. This provides 16 components whose positions are known with milliarcsecond 
precision, thereby creating a testbed for the multipole- Taylor expansion. 

This paper will be organized as follows. The next section describes prior optical and radio 
observations of MG J0414-I-0534, in addition to the latest VLBI map. Section ^ presents the 
formalism for our series expansion, discusses the physical significance of each term, and identifies 
the terms that cannot be constrained due to degeneracies in the lens equation. Some previous 
models for MG J0414-I-0534 are discussed and compared to our technique. Section Q explains the 
numerical methods we employed and presents the best-fit results. Finally, section |5| discusses the 
implications of these results for the mass distribution of the lens and the time delays between 
images. 



2. Observations of MG J0414-h0534 

The radio source MG J0414-I-0534 was first identified by Hewitt et al. (1992) as a gravitational 
lens during a systematic search for lenses in the MIT-Green Bank 5 GHz radio catalog. In both 
optical and radio images it has four components, which have been called Al, A2, B, and C, in 
order of decreasing brightness. The radio components all have spectral indexQ a = —0.80 it 0.02 
( Katz, Moore Sz Hewitt 1997| ), and the optical components are all exceedingly red. However, 
the A1/A2 radio fiux ratio (~ 1.1) and optical fiux ratio 2.5) do not agree, even though 
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gravitational light deflection is independent of wavelength. The discrepancy could be caused by 
a number of factors ( [Angonin-Willaime et al. 1994| ), including dust ( McLeod et al. 1998| ) and 
microlensing (Witt, Mao & Schechter 1995). 



Schechter & Moore (1993) discovered the lensing galaxy in the / band, along with a faint 
object 1" west of component B which they named object X. The optical/infrared spectrum of the 
lensed source resembles a very reddened quasar at Zg = 2.639 =b 0.002 ( p^awrence et al. 1995 ). The 



redshift of the lensing galaxy is zi = 0.9584 it 0.0002 ( p?onry &: Kochanek 1998| ). The best presently 
published optical photometry (uncertainty ~ 0.004 mag) and astrometry (uncertainty ~ 0.02") 
for this system were obtained with Hubble Space Telescope (HST) WFPC2/PC1 observations by 



Falco, Lehar fc Shapiro 1997, who also detected a blue arc extending from Al to A2 to B. 



The deepest radio observations to date did not detect a fifth radio component in the 



system (Katz, Moore & Hewitt 1997). Higher-resolution radio maps, using MERLIN (Garrett 



et al. 1993) and VLBI ( patnaik fc Porcas 1996 ), showed substructure within the four images of 



MG J0414+0534. In particular, the VLBI map resolved images Al, A2, and B into two components 
each, and revealed image C to be extended. Previous attempts to use VLBI constraints to perform 
lens modeling, however, were hampered by the low signal-to-noise ratio of the available data 
( [Ellithorpe 1995| ). 



As part of a 1995 study to test the feasibility of using the NRAO0 Very Long Baseline Array 
(VLBA) to monitor the time-variability of various radio-loud lenses, we obtained 5 GHz images 
of MG J0414-I-0534. These observations and the data reduction procedure are described in detail 



elsewhere ( Trotter 199S ) and will be presented in a future paper; only the results are summarized 
here. The synthesized beam was 1.5 x 3.5 milliarcseconds, and the RMS thermal noise was about 
0.15 mJy/beam, close to the theoretical limit. The peak fluxes in the maps of each of the four 
images ranged from 21 to 110 mJy/beam, allowing a much higher signal-to-noise ratio than was 
previously available. Within each image were detected four components. The compact components 
p, q, and r were labeled in order of decreasing brightness; there is also one extended component, 
labeled s. The four maps are shown in Figure superimposed on a lower-resolution VLA map. 
The sub-components are labeled in the larger maps shown in Figure ^. The locations, fluxes, and 
extents of the components were determined by least-squares fitting to a set of elliptical Gaussian 
parameters. Table |l| lists the locations and fluxes of the 16 VLBI components. 



3. The modified multipole- Taylor expansion 

The motivation behind our method for lens modeling is to expand the lens potential in a 
manner which preserves a measure of mathematical generality and in which the terms that can 
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and cannot be constrained by lensing data are clearly separated. We have tried to arrange for the 
parameters in our expansion to be sensitive to what can be learned from lensing data, rather than 
what one would wish a priori to learn. In this section we will develop this expansion in detail, 
describe the physical significance and constrainability of each term, and compare our expansion to 
other commonly-used parameterizations. 



3.1. The mult ipole- Taylor expansion 

The goal of lens modeling is to deduce the two-dimensional gravitational lens potential, which 
is defined in a way similar to the definitions of Schneider, Ehlers &: Falco (1992) and Narayan & 
Bartelmann (1996): 

m = ^^J dz^M{DLV.z). (1) 

Here r is the angular coordinate measured from an arbitrary optic axis, z is the line-of-sight 
coordinate, and is the Newtonian gravitational potential of the lens. The angular diameter 
distances Dl, Ds, and Dls are from observer to lens, observer to source, and lens to source, 
respectively. As elsewhere in this paper, the speed of light has been set to 1 by a choice of units. 

The lens potential is related to the surface mass density of the lensing object by the 
two-dimensional Poisson equation. 



V?4.(r) = ^.G^air) = (2) 

Implicit in the above equation is the definition of the critical surface mass density, 
Uc = {1/4ttG){Ds / DlDls)- The lens potential determines the image configuration via 
the lens equation, 

s = r - Vr$(r), (3) 

where r is the image position and s is the source position. The typical procedure in lens modeling 
is to adopt a parameterized form for <I>, and then fix the parameters by minimizing an error 
function x^- The error function may represent the deviations between the observed images and 
the images that are projected through the model potential from a model source distribution. 
Alternatively, the error function may represent the deviations between source positions that are 
back-projected from the various images (which is computationally simpler), as explained further 



in section 4.1 



The parameterization we adopt is a multipole- Taylor expansion. The first step is to expand 
the potential $(r) in terms of a complete set of orthogonal basis functions which make the Poisson 
equation separable: 



oo 



^{r,9) = ^o{r) + ^mir) ■ {xcosmO + ysinmO) (4) 

m=l 
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The function ^>o('') is the monopole, and the vector-like functions ^mif) are the higher multipoles. 
Schneider & Weiss (1991) also carried out a multipole expansion of the lens potential in this 
manner. The motivation for this expansion is that each multipole moment of the potential depends 
only on the corresponding multipole moment of the surface mass density. Specifically, the relations 
are (IKochanek 19911): 



rr an(r') an(r') 
$o(r) = constant + 2 In r / dr' r' ^ + 2 dr'r'lnr' ^ ^ (5) 

Jo CTc Jr Cc 

*„(r) = rdr'A^+"^^^^^-— Tdr'r'^^-^^^^^^ (6) 

m Jo ac m Jr cTc 

For smooth mass distributions, only the few lowest-order terms in the multipole expansion 
of the surface mass density are expected to be significant, and the rest may be neglected. By 
reference to equations |5| and ^ it is clear that this is equivalent to neglecting the higher multipole 
components of the potential. For modeling purposes, one may truncate the expansion of the 
potential at the desired level. 

The radial dependence of $o('^) and ^m{r) must also be parameterized to be suitable for 
use in lens modeling. Lensed images constrain V<I> and didj^ at the locations of the images. 
For Einstein rings and quadruple-image lenses ("quads"), the images are typically near the lens's 
characteristic ring radius. Therefore, following the suggestion of Kochanek (1991), we expand the 
radial dependence of each multipole moment of the potential as a Taylor series in the parameter 
p = (r — b)/b, where b is the Einstein ring radius. (The meaning of the Einstein ring radius for a 
non-circular lens is discussed below, in section ^^.) This series should converge quickly for image 
positions located at r ~ 6. For example, \p\ < 0.2 for all of the components in MG J0414-I-0534. 

The resulting expansion of the lens potential is: 

1 1 °° 1 

^Hr,0) = const + + (7) 

- f; ( j + pmM^p + f; i/F„t I • (x cos mO + y sin 

m=l \ I t=2 J 

parametrized by the origin of the expansion {gx,gy), the ring radius b, the monopole parameters 
ft, and the higher multipole parameters M^™, M^^, and F^t- We use t as an integral index, 
contrary to convention, because of the mnemonic value of associating t with Taylor and m with 
multipole. The additive constant /o does not affect light deflection. Because we factor out the 
ring radius b, and then use it as a parameter, we have /i = 1. The m > 1 multipoles with t = 
and t = 1 received the special names M^™ and MJ^^ for reasons that will become apparent in 
the next section. The model parameters gx and gy, which specify the location of the origin of 
coordinates, are implicit in equation |^. To simplify the interpretation of the multipole moments, 
this origin should be centered on the deflector. In section |3.4| it will be shown how this condition 
may be enforced during the model-fitting procedure. 




-7- 



3.2. Physical significance of the expansion parameters 

Although the potential $(r) is the quantity most directly constrained by observations of 
lensing, it is the surface mass density a{r) that is usually of direct astrophysical interest. In this 
section the correspondence between the the parameters in the multipole- Taylor expansion of $(r) 
and the multipoles of cr(r) will be made explicit. This is useful because the multipole moments 
of the surface mass density have a simple physical meaning. This meaning will be reviewed first, 
then the correspondence between the multipole- Taylor parameters and the surface mass density 
multipoles will be given. When a{r) is expanded in a multipole series as in equation it can be 
shown that: 

Mr) = ^rdea{T,e) (8) 



^{r) 



1 r'' 

— / d6a{r,9){'kcosm9 + y smmO). (9) 
TT Jo 

These expressions make the physical meaning of the multipoles clear. The monopole moment 
c7o(r) is the average surface mass density in an infinitesimally narrow annulus of radius r, or, 
equivalently, the angularly-averaged surface mass density. The multipoles <Tm(r) describe the 
distribution of matter around that annulus. In particular, the dipole moment cri(r) points to the 
center of mass Xc of the annulus, viz., 

Xf, 1 (Ti ir) 

— = — 10) 

Likewise, the quadrupole moment <T2(r) arises from an elongated mass distribution. The <T3(r) 
term arises from any triangularity in the mass distribution. Any quadrangularity, e.g. boxiness or 
diskiness, will give rise to a nonzero cr4(r) term. The multipole (Tm{f) can also be expressed as a 
magnitude |<Tm(r)| and an angle V'(t,„('^)i 

(Tm{r) = |cr„(r)|(xcosmV'a™(r) y sinmV'a„(r)), (11) 

rather than by its x and y components. Twisted isodensity contours cause a change in the angle 
i^amif) with radius. In particular, a relative excess of mass at radius r in any of the m directions 
■0o-m(^) + 27rn/m (where n < m is a positive integer) will contribute to the multipole moment 
<Tm(r). It follows from the positivity of the surface mass density that |<Tm(r)| < 2cro(r). The 
multipole moment attains its maximum amplitude only if all the mass in the annulus at radius r 
is clumped at the equally spaced angles tpami^) + ^irn/m, though it need not be evenly divided 
between these angles. In this manner, knowledge of the mass density's multipole moments gives a 
direct picture of the location of the mass. 

We now relate the expansion parameters of the lens potential to the multipole moments of 
the surface mass density, by expanding equations || and ^ in a Taylor series about r = b. For the 
monopole, 

- ^^*|?^o(r) , (12) 
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of which the first four terms are: 



r=b 



r=b 



r=b 



const 



-1 + 



2ao(6) 



CTc 



+ 2 6 



d_ f cro{r) 
dr 



(13) 
(14) 

(15) 
(16) 



r=b 



As mentioned previously, the constant term /q can be ignored because a constant offset in the 
potential does not affect light deflection, nor does it affect the differential time delay between two 
images of the same source. (Although /o does affect the propagation time of the light from source 
to observer, only the relative or differential delay between two different images can be measured.) 
The linear term, /i, sets the ring radius h and does not have any other adjustable coefficient. In 
fact, equation |l^ defines the Einstein ring radius h not only for this expansion, but also for a 
general potential. It can be shown that the average surface density within the radius b so defined 
is the critical surface density cic ( Schneider, Ehlers &: Falco 1992| ). The quadratic term /2 gives the 
angularly-averaged surface mass density at the ring radius. In the next section we will show that 
this term cannot be constrained by lensing phenomena, due to the mass-sheet degeneracy. The 
cubic term is therefore the first constrainable term yielding information about the radial profile 
of the mass distribution. 

The correspondence between the higher multipoles of '&(r) and (T(r) are obtained by 
differentiation of equation ^ Before presenting the first few terms of this correspondence, we 
define two quantities and B^, which are attributable to the m**^ multipoles of the mass that 
is, respectively, exterior and interior to the Einstein ring radius: 



A, 



m Jb Uc 



m 



dr r 







(17) 
(18) 



With these definitions, the correspondence between the higher multipoles may be written 
compactly. In general, for m > 1, 



mt 



1 f/* 



r=b 



Explicitly, the first four terms of the Taylor expansion are: 

1 



M!:^"' = F 



mO 



62 



^m{b) 



Am + B^ 



(19) 



(20) 
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wM^'*^ = F 1 = —h—^ 

mivi^ — ' ml — ip. dr 



1 , d 



(r) 



1 



F„2 = - ^h'^^m{r) 



^ - ^b^^^mir) 



r=b 



m{Am - Bm) 



m{m — 1)A^ + m{m + 1)B^ 



(21) 
(22) 



m{m — l)(m - 
I 2o-^(6) 



2)A^ - m(m + l)(m + 2)B^ 

d CTm{r) 



2 6 



CTc. 



(23) 



Since the parameter of the t = term is the sum of the exterior and interior multipoles and 
B^, we label it M^™ rather than F^o- Likewise, the parameter of the t = 1 term, labeled M^^, 

is the difference between the exterior and interior multipoles. Higher-order terms {F„it for t > 2) 
depend on A^^ for t ttt-, cind on B^j^^ for all t, as well as on the behavior of crj^i(^T^ near r — b. In 
the next section it will be shown how to modify this expansion in order to separate explicitly the 
effects of and B^ from the contributions that depend on <Tm ('") ill the vicinity of r — b. 



3.3. Separation of internal and external contributions 

The physical meaning of our parameterization becomes clearer if A^ and B^ are used directly 
as parameters, instead of M^™ and M^*'^. Furthermore, the sum over t of all the terms involving 

and B^ can be calculated exactly, leaving only the effect of the mass near the ring radius left 
in the Taylor series. When this is done, the resulting "modified" multipole-Taylor expansion is: 



1 1 °° 1 

p$(r,0) = const + p+-p^f2 + Y,-^p' ft 



t=3 



+ E 

m=l 



(24) 



Bm + X] Jiy^rnt \ ' (x COS 777.6' + y sin 7776 



t=2 



In this expression, the terms G^t are the portions of F^t that are left when the contributions 
proportional to A^ and B^ are subtracted out. The generic form for Gmt is {m > 1, t > 2): 



1 J* 

G + — --?)*_*'i«'arrmg/ > 



r=b 



where 



^nearring/ X ^ - I - 

^' 777 U 



1 /r\" r'' fr'V"^ dr' amir) 1 



r 

777 \b 



/X 1+m 



dr' (Tmir) 



Or. 



(25) 



(26) 



Explicitly, the first four terms in the Taylor expansion are: 



Gmi = 



(27) 
(28) 
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G 



m2 



'm3 



rib) 



0"c 

d CTm{r) 



(29) 

2b^^^^^ -2:1^ (30) 

dr CJc r=b CTc 

We see that these parameters depend only on the distribution of the deflector mass near the ring 
radius. In particular, Gm2 gives the m^^ multipole moment of the matter located at the ring 
radius itself. 



3.4. Degeneracies of the lens equation 

Certain parameters in the multipole- Taylor expansion cannot be constrained by observations 
of gravitational lensing, due to the so-called "degeneracies" of the lens equation. A degeneracy, as 
explained by Gorenstein, Falco &: Shapiro (1988), is a family of distinct lens potentials and source 
distributions that all produce the same image configuration. Thus the observation of a given image 
configuration does not permit the actual lens potential and source distribution to be deduced, but 
rather only the degenerate family to which they belong. Breaking the degeneracy requires either 
additional assumptions or direct observations of the deflector itself (such as velocity dispersion 
measurements; see Falco, Gorenstein &: Shapiro 1985| ). 



The unconstrainable parameters must be fixed during the minimization of the error function 
X^, or else convergence is impossible. This section will identify those parameters and the fixed 
values we have chosen for them. The first degeneracy involves the choice of origin. From 



equations 18 and U follows, for m = 1, 



B, = U\rr^^ = ^, (31) 

0'^ Jo CTc b 

where Xcom is the center of mass of all the mass interior to b. This parameter changes with 
the choice of origin, but obviously the arbitrary choice of origin cannot affect the image 
configuration. Therefore one may allow either Bi or the location of the origin to vary, but not 
both simultaneously. We decided to impose Bi = and allow the location of the origin to vary 
during the minimization procedure. This choice forces the center of the interior (r < b) mass to 
be at the origin, thereby simplifying the physical interpretation of the higher multipoles. 

The second degeneracy is the "prismatic" degeneracy: any term of the form c • r in the 
potential, where c is a constant vector, is unconstrainable ( |Falco, Gorenstein &: Shapiro 19"85 



Gorenstein, Falco fc Shapiro 198^ ). The contribution of the exterior dipole Ai to the potential 



(equation |^ or equation 24) is the additive term — 6A1 • r, which is precisely a prismatic-degenerate 
term. Consequently, we may arbitrarily set Ai to zero during the model-fitting procedure. (If the 
expansion in equation |^ is used, rather than the "modified" expansion with internal and external 
contributions separated out, then the center-of-mass degeneracy and the prismatic degeneracy can 
be taken into account by omitting the M|"™ and M]*'^ terms, which is equivalent to setting both 
Ai and Bi to zero.) 
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The last and most problematic degeneracy is the "mass-sheet degeneracy" ( [Falco, Gorenstem 



fc Shapiro 1985| ; Gorenstein, Falco &: Shapiro 1988 ). Two potentials $(r) and $'(r) cannot be 



distinguished by the observation of lensed images if they are related by the transformation: 

$'(r) = (1 - K)Hr) + ^. (32) 

This is equivalent to reducing the surface mass density by a factor (1 — k) and adding a sheet of 
uniform mass density nac- 

cj'(r) = (1 — k) cr(r) + nac- (33) 

The corresponding relations between the multipole- Taylor parameters of the original and 
transformed potential are: 



t > 3 



/o = (1 - ^)/o + \'^ = ^o'^^*' 


(34) 


/( = /l = 1 


(35) 


= (1 - k)/2 + K 


(36) 


ft = (1 - «)/t 


(37) 




(38) 




(39) 


= (1 - K)Bm 


(40) 


G'mt = (1 ~ K)Gmt 


(41) 



The constant offset has been altered, but this is immaterial since adding a constant to the potential 
has no effect on light deflection or differential time delays. The transformation must leave the 
ring radius unchanged, so the coefficient of p in the monopole part of the potential is unchanged. 
The quadratic monopole term, /2, which depends on the surface mass density at the ring radius 



(Eq. |15D has been changed, because the surface mass density at the ring radius has been changed. 
All the coefficients of the other terms in the model have been scaled by a factor 1 — k. These 
simultaneous adjustments of all the ft (t > 2), Am, B^, and Gmt parameters leave the model 
predictions unchanged. Therefore one of these parameters should be fixed in some manner during 
the model-fitting procedure to permit convergence of the minimization algorithm. 

We note that the transformation of equation ^ does affect the time delays, 

At' = (1 - K)At. (42) 

This however does not permit observations of the time delays to be used to determine the scaling 
factor (1 — k), since a rescaling of the time delay will be confused with a rescaling of the Hubble 



parameter. We discuss this issue further in section 5.2 



To break the degeneracy, and permit convergence of the minimization algorithm, we may 
arbitrarily select a value of f2', by equation |l^, this is equivalent to an arbitrary choice of o"o(6). 
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We find it convenient to fix = 0, which happens to be correct for an isothermal sphere, for 
which (To(6) = fTc/2. If this is not true of the actual mass distribution, then the fitted (primed) 
amplitudes of all the m ^ multipoles and the t > 3 monopole parameters will have been scaled 
by the same factor 1 — k relative to the true (unprimed) amplitudes, where 



1 



2 1 



0"c 



(43) 



These scaling factors are tabulated in Table ^ for several different choices of the actual monopole 
potential. For example, suppose the choice /2 = is made, but the true potential is actually a 
point mass rather than an isothermal sphere. Then, in reality, $o(r) = 6^1nr, and o"o(6) = 0. 
All of the best-fit parameters (//, A^, etc.) and the model time delays (At') are smaller than 
those {ft, Am, At, etc.) of the actual mass distribution, having been scaled by the same factor 
(1 — k) = 1/2. Note that this rescaling only affects the magnitudes of the multipoles; it does not 
affect their angles. It also does not affect the center-of-mass parameters gx and gy, or the ring 
radius b. 

Once we have set /2 = 0, the next term in the monopole expansion, /g, is the lowest-order 
monopole parameter available to give information on the radial distribution of mass. After 
rescaling, its amplitude is related to the true surface mass density by 



4 In 

dr 



1 



r=b 



(44) 



in contrast to equation 16. It depends on the fall-off of the surface mass density near the Einstein 
ring radius, and is most sensitive to that falloff when o"o(6) is close to cJc- For surface mass densities 
that do not increase with radius, < 1. The values of for various monopole potentials are 
given in Table |5[ The lowest value for /g listed in the table is —1, for a mass sheet, but can 
be even more negative for potentials with a core radius, or if o"o(r) drops abruptly near the ring 
radius. 

The mass-sheet degeneracy also affects the interpretation of the amplitude of the quadratic 
multipole parameter G^2) whose fitted value is related to the true surface mass density by 



m2 



Mb)' 



(45) 



in contrast to equation ^ Non-negativity of the surface mass density limits the amplitude of 

g: 



m2- 



2Mb) 



CJc - ao[b) 

This limit ranges from zero, for cto(&) = 0, to +oo, for croib) = 0"c. For a mass distribution with a 
singular isothermal sphere monopole, the limit is |G^2l — 2., For mass distributions more centrally 
concentrated than a singular isothermal sphere, the limit on |G^2l would be lower. The limit is 
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only attained if all the mass at the ring radius is clumped into point-like perturbers in any of the 
m allowed directions. 

In summary, after expanding the lens potential in a multipole- Taylor series, then explicitly 
separating the contributions from the internal and external multipoles, and then fixing several 
parameters to zero (as described above) because of the degeneracies in the lens equation, the final 
parameterized form of the potential is: 

1 °° 1 



^ {(1 + prA'^ + (1 + pr'^B'^} ■ (i cos me + y sin 



mb 



m=2 
oo oo 



.1 P ^mt ' ("^ ^os mO + y sin m6 

m=l t=2 



(47) 



Here and elsewhere, whenever the distinction is important, we have used primes to identify 
fitted model parameters, reserving the non-primed symbols for the parameters that describe 
the actual gravitational potential. These two sets of parameters differ because of the arbitrary 



degeneracy-breaking choices described in this section. Equation |47| is the parameterization we 
used to compute best-fit models for MG J0414-I-0534 (although, in a few cases, we tried using 
the parameter M™™' from the original multipole- Taylor series, equation |^ in place of and 
B^). A summary of the physical significance of each term in this expansion is given in Table ^. 
Since the analytic forms of the first and second derivatives of the potential are used during the 
model-fitting procedure, these derivatives are presented in Tables ^ and ^. 



3.5. Comparison to other parameterizations 

Any parameterized form of the lens potential <&(r) can be compared to ours by expanding it 
in a modified multipole- Taylor series. In this section we carry out this procedure for a few forms 
for the lens potential that are widely used. 

Often the lens potential is taken to be a combination of a monopole term and a quadrupole 
term, just as Kochanek (1991) did. The particular form of each of these is usually prescribed 
by either physical preconceptions or ease of computation. In Table |5| we compare five possible 
monopole potentials. These all give rise to the same linear term which sets the Einstein ring radius. 
The first significant term that distinguishes them is therefore f^, since /2 is not constrainable 
due to the mass-sheet degeneracy. If the data are not sufficient to distinguish among the various 
possibilities for /g, then there is little point in using more complicated forms for the monopole 



term, such as de Vaucouleurs radial profiles (Ellithorpe 1995) 



There are three different quadrupole terms that are commonly used to accompany the 
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monopole term. The first is an external quadrupole, of the form: 

c^Q(r,0) = V cos 2(0-0^) (48) 

This term is sometimes referred to as an "external shear." In our scheme, this is equivalent to 
choosing a non-zero value of A2, with magnitude IA2I =7/2 and angle ipA2 = vr/2 + 6*^. (Our 
notation has ipA2 '^A2 + as the angles to the mass excess, whereas 9-y and 9-y + tt are the 
angles to the mass deficit.) The second is an internal quadrupole term, parameterized as: 

^Q{r,9) = -^^cos2{e-e^), (49) 

which is equivalent, in our scheme, to choosing a non-zero B2 with magnitude [B2I = 7/2 and 
angle ipB2 = 0^. Finally, there is the case of a "mixed" quadrupole, 

^Q{r,e) = -^rcos2{e-e^), (50) 

which is obtained by truncating the potential of a singular isothermal elliptical potential at the 
quadrupole term. To order p, this is equivalent to the choice IA2I = 87/8, IB2I = 7/8, and 

= "0-82 = ^7- It is worth noting here that A2 + B2 causes tangential image displacements for 
images near the ring radius, whereas A2 — B2 causes radial image displacements. The balance 
between the internal and external portions of the quadrupole therefore determines the radial 
displacement that accompanies the tangential displacements caused by the quadrupole moment of 
the mass distribution. The choice of a "mixed" quadrupole is essentially a particular choice for 
this ratio. 



3.6. Effects of external perturbing masses 

If there is a deflecting object along the line of sight to the source besides the primary 
lensing galaxy, this secondary deflector can be modeled in one of two ways. It could be treated 
separately from the primary deflector, with extra parameters for its location and mass distribution. 



Alternatively, the parameterization of equation ^ could be used alone, so that the influence of the 



secondary perturber would be reflected in the values of the multipole moments. Since the location 
(or even the presence) of a perturbing object is not known a priori, it is useful to compute the 
effect on the parameters of the multipole- Taylor expansion that would be caused by a perturbing 
object far removed from the primary deflector. 

A perturbing point mass located further from the origin than any of the lens images 
contributes to the external multipole moments to all orders m, with amplitudes 

I fh\'^ fhp\ 
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where b is the ring radius of the principal deflector, 6p is the ring radius of the perturber, and R 
is the distance to the perturber. This can be derived by expanding its delta-function surface-mass 
distribution as a multipole- Taylor series. The m = 1 term has no effect, due to the prismatic 
degeneracy, so the dominant term is the quadrupole moment, which is often called the external 
shear. For small b/R, it should be adequate to represent the effect of the perturber by the first 
few multipoles. Conversely, the ratios of the fitted amplitudes (if they can be attributed solely 
to a perturbing influence) permit the distance to the perturber, and its strength, to be deduced. 

How do radial and tangential extent in an external perturber affect the external multipole 
amplitudes? Consider a perturber with its center of mass located at a distance R from the center 
of the principal deflector, with surface mass density uniform over a region of radial extent AR 
and tangential extent RAO. Such a perturber contributes to the external multipole moments of 
the potential to all orders m in the multipole expansion, with amplitudes 

1 fb\"' fbp\^j ^ , m(m + l) fAR\^ m{m + 1) 
'^'"' = mUJ [t) y ^ 24 [-If) 2^(^^) 

+ O (^)' + O (A^)^ + O {AOf O (^)' } , 

(52) 

The net effect, for perturbers with similar radial and tangential extents {AR = RAO) is that 
there is no effect through O (■^^ • Therefore, using the formulae for a point mass perturbation 
(eq. |5l|) should cause no problem for a perturber that is located far away as compared to its 
extent. Unfortunately, for a nearby extended perturber, the interpretation of the fitted external 
multipole amplitudes is less simple. 



4. Application to MG J0414-|-0534 

A model of a gravitational lens consists of two parts: a model of the surface brightness of the 
source, as it would appear in the absence of lensing, and a model of the lens potential. An ideal 
model reproduces the observed image of the system, pixel for pixel, by mapping the source surface 
brightness through the lens potential and then convolving with the detector response. This is the 



aim of modeling techniques such as LensClean ( Kochanek &: Narayan 1992 ; Ellithorpe, Kochanel: 



fc Hewitt 1996 ), which use the information from every pixel to constrain the model. 



It is computationally faster to employ a much smaller number of constraints that capture the 
most important features of the observed image. Multiple images of nearly-pointlike sources, such 
as the 16 components seen in MG J0414+0534, are succinctly described by the locations and flux 
densities of the pointlike components. With this simplification, called "point modeling," a much 
wider range of models can be explored. 



We chose to use only the positions of the 16 components of MG J0414+0534, and not their flux 
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densities, as model constraints. There were two reasons for this. One, the relative uncertainties in 
the fluxes are much larger than the relative uncertainties in the centroid positions. Therefore the 
inclusion of flux information would not contribute much to the error statistic ^ . Two, there are 
significant discrepancies between the radio measurements and optical measurements of the flux 
density ratios of components in MG J0414+0534. This discrepancy has been variously attributed 
to microlensing ( Witt, Mao &: Schechter 1995| ), variable extinction ( McLeod et al. 1998 ), and/or 
substructure in the lensing galaxy ( Mao fc Schneider 1998 ). Some of these effects may not be 
significant for radio observations, but we chose not to use flux density information at all. 



4.1. Model- fitting algorithm 

Our constraints consist of the observed positions r^fj, where the index a runs over the 4 
different images (Al, A2, B, and C) and the index k runs over the 4 components (p, q, r, and s) of 
each image. The model consists of presumed source positions for each component, along with 
the mapping rQ;(s/[;) provided by the lens equations (eq. |3]) using our expansion of the potential 
(eq. |4^) truncated to a prescribed order. 

Assuming that the observational errors are Gaussian, the maximum-likelihood values of the 
parameters can be determined by minimizing the familiar chi-squared statistic, 

4 4 

= (ra(sfc) - rfca) • S^^^ • (r„(sfc) - rfe„), (53) 

k=\a=\ 

where S^q is the error covariance matrix. Because the lens mapping is much easier to apply in the 
direction r — > s than the reverse direction, it is of great computational advantage to compute an 
approximate value of by evaluating the model errors in the source plane rather than the image 
plane ( |Kayser et al. 1989 ). By Taylor expansion, the displacement between the model source 



position and the actual source position is 

Sfc - &(vka) = M"^(rfca)(ra(sA;) - Ffc^) + {(Va{'&k) " Tfca) • VrM"^ (rfc^)) (rQ,(Sfc) - rfc„) + . . . , (54) 

where M~-'-(rfcQ) is the inverse magnification matrix. 

For good models, the difference between the modeled and observed image-plane positions, 
ra(sfc) — rfcQ, is small, so the higher order terms in the expansion may be neglected. Equivalently, 
the change in magnification between the observed image location and the model-predicted image 
location is assumed to be negligible. The resulting "source-plane" approximation for is: 

4 4 

~ = E E (^fc - s(rfca))M(rfc,) • S^^^ • M(rfc„)(sfc - s(rfc„)) (56) 

k=\a=\ 
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This approximation is useful because it is valid near the chi-squared minimum, and we are 
unconcerned with the behavior of the function far from the minimum. The global minimum of 
both and Xs is zero, which occurs when the model exactly reproduces the observation. We 
expect that even with noise and measurement error, Xs has a global minimum corresponding to 
that of the true chi-squared, and that no lower minimum is introduced by this approximation. 
This approximation fails for very poor deflector models that do not reproduce the image locations 
at all — but this is of little concern, since the high would cause the model to be rejected 
anyways. The approximation can also fail if the error in the observed image locations is large 
enough to encompass a region in which the magnification matrix varies significantly. Because of 
this danger, after the minimum of was found for each model, we computed the true to check 
the source-plane approximation. For models that adequately satisfied the observational constraints 
(as described below), typically differed from x^ by only 0.2 — 0.6%, and at most 1.5%. Since 
this is much less that the x^ increment used to find confidence limits on model parameters, the 
source-plane approximation introduced no appreciable error. 

Since the first and second derivatives of our parameterized potential are available in analytic 
form, the magnification matrix and thus Xs are easy to compute. These derivatives are listed 
in Tables |3| and Furthermore, since xi is quadratic in the source positions s^, the optimal 
source positions are easily computed for given values of the model parameters. Consequently, 
the numerical minimization need not include a search through the source positions in addition to 
the parameters of the deflector model. This reduction in the number of dimensions of the search 
space permits a vast computational speed-up. (If fluxes also are used as model constraints, a 
corresponding source-place approximation can still be made. The resulting is quadratic in both 
model source positions and model fluxes, so the optimal model positions and fluxes can still be 
found analytically Trotter 1998| .) 

To perform the minimization of Xs we employed a variant of simulated annealing described 
in Press et al. 1992. Simpler methods, such as a straightforward downhill walk via the Powell 
direction set method, or the downhill simplex method ( press et al. 199^ ), encountered difficulties 
with local minima. The initial value for the center of mass was set to the galaxy position 
observed in the HST images of Falco, Lehar & Shapiro (1997). The initial value of the ring 
radius b was taken from the previous best-fit models of Ellithorpe (1995), and the initial values 
of all m > 1 multipoles were set to zero. To explore the region of parameter space near these 
physically- motivated starting values, we used a moderately low starting temperature equal to 
1% of the initial value of Xs- High starting temperatures allow the possibility of escaping the 
physically reasonable portion of parameter space and becoming trapped in deep and distant local 
minima. 

The higher multipole parameters were represented by their Cartesian components, e.g. 
X • and y • Ag, rather than by amplitude and angle. The potential depends linearly on 
these parameters, and thus the chi-squared depends roughly quadratically on them (far from the 
minimum, at least), allowing for a more robust minimization. However, confidence limits were 
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computed for the amplitude and angle of each multipole parameter, rather than its Cartesian 
components, because the amplitude-angle representation is more useful for visualizing the mass 
distribution, and because amplitudes are affected by the mass-sheet degeneracy but angles are not. 

For the value of to be used to calculate confidence intervals for the model parameters, 
the estimates of the observational errors (as represented in the covariance matrix Ska) must 
be accurate. Unfortunately, it is difficult to make accurate error estimates of VLBI centroid 
positions, because of the complicated and nonlinear process of deconvolution by "cleaning" and 
self-calibration. In addition, if there are magnification gradients across the image, the image 
centroid may not be exactly the image of the source centroid, even though the point-modeling 
approach assumes so. The separation between the image of the source-centroid and the centroid 
of the image is approximately the angular size of the image multiplied by the fractional change 
in the magnification over the extent of the image. (The expression for the discrepancy is given 
by Trotter (1998), along with a correction to Xs to account for it. We did not use this correction 
because it requires the calculation of the third derivatives of the potential and accurate estimates 
of the intrinsic source size.) 

For these reasons we report the results using three different methods to estimate the positional 
error in each component. The first, a crude upper limit, is the image size convolved with the 
VLBA beamwidth, which we call "fit-size errors." The second estimate for the positional error is 
a lower limit: the statistical error in the centroid position due to thermal noise in the map. We 
report this as "statistical error." The third estimate is the quadrature sum of the statistical error 
and the width of the deconvolved image (the intrinsic component size). This "stat-width error" 
makes some allowance for magnification gradients as well as deconvolution error. We believe the 
stat-width estimate to be the most accurate of the three estimates, but since this judgment is not 
rigorous we report the results for using all three estimates. Confidence limits on each model 
parameter were determined using the "stat-width error," by stepping the parameter away from 
the Xs minimum, while minimizing over all other parameters, until the Ax'^ appropriate for 68.3% 
confidence limits was obtained. 

4.2. Model results 

Table ^ summarizes the goodness-of-fit for a variety of model potentials. For each model the 
number of parameters and number of degrees of freedom are listed, along with the minimum x^ 
obtained using each of the three different estimates of positional errors. In this section we review 
these results. 

Before using the modified multipole- Taylor expansion, we tested three simpler 
parameterizations that have been used in previous attempts to model MG J0414-I-0534 
based on lower-resolution radio and optical data. These simpler potentials were a singular 
isothermal sphere plus external shear (SIS-I-XS), a point mass with external shear (PM-I-XS), and 



- 19 - 



a singular isothermal elliptical potential truncated at the quadrupole moment (STEP) (see e.g. 
Hewitt et al. 1992| , [Falco, Lehar fc ShapiroT997| , [Ellithorpe 1995| ): 



^'sis+xs(r,^) = hr + ]pr^ cos2{e - e^) (57) 
$PM+xs(r,^) = 6^ In r + ^7r2 cos 2(0-0^) (58) 



^siEp(r,^) = 6r-^r cos 2(^-0^) (59) 



The results for these 5-parameter models are shown in the top three lines of Table ^. They 
are all very poor fits. Qualitatively they fail to reproduce the detailed VLBI structure of the four 
components. Quantitatively, even for the upper-limit (fit-size) errors, the is more that 3 x 10^ 
standard deviations away from the value that would be expected for observations matching these 
models. 

The rest of the models in Table ^ are multipole- Taylor expansions truncated in various 
ways. They are labeled with symbols indicating the terms that are present in the expansion. 
The dominant terms in the multipole- Taylor expansion — the only terms that cause shifts in 
the image positions for images at the ring radius — are the first (linear in p) monopole term, 
6, and the first two (constant and linear in /o, or external and internal) m > 2 multipole terms, 
Am and B^. Accordingly, all of the models included the b term, which sets the Einstein ring 
radius, as well as the internal quadrupole (B2) and external shear (A2) terms which account for 
ellipticity of the deflector mass distribution and the dominant effects of any external perturber. 
The model containing these three terms and no others is given the schematic name b + A2 -|- B2 , 
and is the fourth model listed in Table ^ It differs from the SIS-I-XS model only by the addition 
of the internal quadrupole B2 term. The addition of this second quadrupole term causes a vast 
improvement in the fit; the minimum is lowered by two orders of magnitude, although the 
model is still not in formal agreement with the data. 

All of the other models are labeled with a sign and the terms they contain in addition to 
the three terms 6, A2, and B2. The next batch of models, as indicated in Table ^, each include 
only one term in addition to these three. The next most significant terms in each of the m = 
through m = 4 multipoles were tried. The m = 3 multipoles were considered because they are 
the next terms in which effects of an external perturber would appear, and would also account 
for any lopsidedness in the mass distribution of the lens galaxy. The m = 4 multipoles were 
considered to account for diskiness or boxiness of the lens galaxy. Higher multipoles, m > 5, were 
not tried, as there was no physical reason to expect them to be significant. Also tried were the 
mixed-internal-and-external terms, M™™ and M™™, of the original multipole- Taylor expansion, 
equation ^. Using M™™ instead of or B^ strikes a different balance between the internal and 
external contributions to the m}^ multipole. 



In all cases the fits were an improvement over the three-term model, but the best results were 
obtained by adding either an external octupole A3, or a mixed-internal-and-external octupole 
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Mg"™. These two models fit the image positions weh enough to satisfy the fit-size (upper limit) 
position errors, though not well enough to satisfy the two tighter error estimates. For both of 
these models, the angle of the external quadrupole (74.4° ± 0.2° E of N for +A3; 74.1° ± 0.2° for 
-I-M3"™') is consistent with the optical isophote angle of the deflector as observed in the WFPC2 



image of Falco, Lehar & Shapiro 1997 (71° it 5°). For the -I-A3 model, the direction of the internal 
quadrupole (75.4° ± 0.5° E of N) also agrees with the WFPC2 optical isophotes. The best-fit 
parameters of this 9-parameter model (-I-A3) are displayed pictorially in Figure ^ The -t-MI*^™ 
model has a somewhat smaller value of although in this case the internal quadrupole B2 
(87.4° ± 0.9° E of N) and the isophotes are misahgned by 16°. 

To each of these promising models, +A3 and -l-MI*^™, was added the next radial term in each 
m < 4 multipole, one at a time, as indicated in the third block of Table ^. These models had 10 
or 11 parameters, depending on whether the additional term was monopole or not. For the sake 
of comparison, a model with 11 parameters but employing m = 4 multipoles instead of m = 3 
multipoles was also tried (and fared very poorly). 

The models including the external octupole A3 outperformed the models including the mixed 
octupole Mg^™. The model with the lowest was +A3 -|- A4 (external m = 3 and 4 multipoles), 
followed by the model +A3 -|- B3 (external and internal m = 3 multipoles) Both of these models 
oversatisfy even the stat-width error estimates (indicating, perhaps, that these error estimates may 
be too large), though neither model is formally consistent with the lower-limit error estimates. The 
best-fit parameters for the most successful 11-parameter model, +A3 + A4, are shown pictorially 
in Figure ^. The implications of the success of this and other models will be considered in the 
next section. 

The lowest-order term in the multipole- Taylor expansion that is both constrainable and 
sensitive to the radial distribution of mass is f-^. In the next round of modeling we added the 
term /g to each of the models of the previous group, in order to determine whether the radial 
dependence of the potential could be usefully constrained. The results are shown in the last 
block of Table ^. The fitted parameter values have large error ranges (e.g. = 1.131q52 for the 
model +A3 + A4 + /3) compared to the range of interesting values, which extends from f^ = 
(singular isothermal sphere) to /3 = 1 (point mass). More problematic is that the value of this 
parameter depends sensitively on the presence or absence of the other multipole components 
in the model, with values ranging from —9.7 to 2.3 for 12-parameter models that adequately 
satisfy the stat-width error estimates. It is clear that useful information on the radial profile of 
MG J0414+0534 is unavailable from this data. 



^ The model +A3 + G22 produced the unphysical value IG22I = 0.5. This would require that at least 25% of the 
mass in a narrow annulus at r = 6 be involved in driving the quadrupole term — or more, if mass is not concentrated 
into points. (This assumes that the mass of the deflector is not more extended than a singular isothermal sphere.) See 
equation |46| To produce sizable image displacements, IG22I must be large because it must overcome the small factor 
p. This term is apparently compensating for the internal and external quadrupoles, which are displaced relative to 
their orientations in the +A3 model. 
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5. Discussion 

5.1. Implications for the mass distribution 

The best-fit model of the lens potential with 11 parameters included an m = 3 external 
multipole and an m = 4 external multipole. What are the implications of this success for the mass 
distribution of the deflector? To recapitulate, the explicit form of the model potential in this case 
is: 

= b'' (^p- (^{l + pfA'2+ B^^ • (x cos 29 + y sin 29) - (60) 

(1 + p)^ A3 • (x cos 39 + y sin 36*) - (1 + p)^ A4 • (x cos 46 + y sin 46*) 

where 6 is measured north of west about the center of mass {gx^Qy) as given by 
ianO = {y — gy)/{x — gx), and the radial parameter p is given by 

P={^^{^-9.y + {y-9y?-b^ /b. (61) 

Tabled contains a list of the best-fit parameters. In addition, Table |8| lists the image magnifications 
predicted by this best-fit model (which may be compared to the flux ratios in Table ||). Figure ^ 
shows both the observed and modeled image locations for each of the sixteen components of the 
VLBI map, along with the stat-width error ellipses. 

In this model, the center of mass of the mass distribution interior to r = & is located at 
Aa = -1".0788 ± 0".0020, A6 = 0".6635 ± 0".0012 relative to the correlation center at component 
Alp. Since the position of component Cp relative to the correlation center was Aa = — 1".9454, 
AS = 0".3004 (with position errors negligible compared to those of the model galaxy, see Table |^) 
the model center of mass is Aa = 0".8666 ± 0".0020, A6 = 0".3631 ± 0".0012 relative to Cp. This 
is in agreement with the optical centroid of the lens galaxy as observed by Falco, Lehar &: Shapiro 
(1997) even though the optical position was not used as a modeling constraint. The observed lens 
galaxy position in i?-band, relative to component C, was Aa = 0".90 ± 0".05, A6 = 0" .32 ± 0".05; 
the observed position in /-band was Aa = 0".86 ± 0".05, AS = 0".36 ± 0".05. 

As is apparent in Table ^ and Figure ^, the directions of mass excess implied by both the 
internal and external quadrupoles (Ag and Bg) agree with the direction of the optical isophotes 
(71° ± 5°) observed by Falco, Lear Sz Shapiro (1997). Interestingly, one of the directions of mass 
deficit implied by the external m = 4 multipole, A4, is also in agreement with the optical isophote 
angle. The eastern direction of the external octupole, A3, is aligned with the external quadrupole 
moment within 6°, and is within 10° of the observed isophote angle, although in neither case do 
the formal confidence regions overlap. It is possible that the alignments between the multipole 
angles and the optical isophotes are coincidences. If the isophote angle were selected at random, 
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the chance of agreement with the direction of mass-excess indicated by would be 6%, given the 
quoted confidence ranges. The chance of agreement with the direction of either the mass excess or 
mass deficit implied by A4 would be 27%. The chance that A3 and the isophote angles would be 
as closely aligned as they are is 34%. 

Bearing this in mind, we entertain three speculations regarding the origin of the external 
multipoles in the best-fit model b + A2 + B2 -|- A3 -|- A4: 



1. All of the external multipoles are attributable to the mass distribution of the lens galaxy. 
This explains the alignments of the various multipole angles with the optical isophotes, but 
it implies that no external perturber (i.e. neither object X nor the group of galaxies to the 
southwest) contributes significantly to the features of the potential we have modeled. A 
singular isothermal elliptical potential would have a ratio of external to internal quadrupole 
of IA2I/IB2I = 3, which is consistent with the value 2.90 it 0.17 obtained for this model. The 
ellipticity (1 — b/a) of the isopotential contours near the ring radius for the fitted model 
is e$ = 2IA2 -|- B2I = 0.120 lb 0.002. In contrast the ellipticity of a singular isothermal 
elliptical potential (SIEP) having an isodensity ellipticity of = 0.20 it 0.02 (equal to the 



mean ellipticity of the fitted isophotes of [Falco, Lehar fc Shapiro 19971) would be e$ = 0.07. 
However, the non-zero A3 implies that the outer galactic halo is asymmetric, with more mass 
concentrated near one end of the isophote axis than the other. Furthermore, the magnitude 
of A4 (4.1 X 10~^) is an order of magnitude larger than the value it would assume for an 
singular isothermal elliptical potential (6 x 10~^). The direction of A4 implies that the mass 
distribution is box-like, rather than disk-like as for a SIEP. 

2. A2 and A3 indicate an external perturbing mass to the east. In this case the alignments 



of all multipoles with the isophotes are accidental. According to equation 51, a point mass 
with an Einstein radius of 0.95 ± 0.05 arcseconds, located 3.2 it 0.2 arcseconds away, would 
supply approximately the appropriate values of A2 and A3, but would not account for A4. 



No perturber of any kind is seen to the east in the optical images of Falco, Lehar & Shapiro 



1997 



3. A3 and A4 indicate an external perturbing mass at about —65° it 4° N of W. In this case, 
as above, the alignments of all multipoles with the isophotes are accidental. A point mass 
2.24 lb 0.24 arcseconds away, with Einstein radius 0.56 it 0.10 arcseconds, would supply the 
proper IA3I and IA4I, but would also make a significant contribution to A2. The residual 
component of A'2, which would presumably be due to the lensing galaxy, has magnitude 
0.06 and direction 1.5° or 181.5° north of west. This is somewhat problematic because the 
residual quadrupole does not agree in angle with the optical isophotes, and its magnitude is 
uncomfortably large (IA2I/IB2I = 3.8). However, there is an object seen in the optical image 



of Falco, Lehar & Shapiro 1997 about 5" away in the direction —54° N of W. 



None of these speculations is entirely satisfactory. The first speculation attributes a peculiar 
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shape to the galactic halo; the second invokes an external perturber that does not seem to be 
present; the third compels the galaxy to produce an unusual quadrupole moment. We favor 
the first interpretation, because it is hard to arrange for an external perturbation to produce 
A3 without ruining the suggestive alignment of A2 with the isophotes, but admit that this 
interpretation is debatable. 

This ambiguity of interpretation illustrates both the appeal and the frustration of the 
multipole- Taylor technique for modeling lens potentials. The technique makes few preconceptions 
about the shape of the potential, which in principle may lead to unanticipated discoveries about 
the mass distribution of the deflector, including the shape of the halo (of which little is presently 
known). However, precisely because of this mathematical generality, there is no determinative 
way to correlate features of the potential with observed astrophysical objects (the lens galaxy or 
external perturbers). 

Finally, we comment on the unconstrainability of the radial distribution of the mass monopole. 
It is not terribly surprising that we were unable to usefully constrain the parameter /g, the 
lowest-order parameter containing information about the radial distribution. A large change in 
/g is needed to cause a small change in the radial positions of images near the ring radius. By 
contrast, even small values of A^ or (for m >2) cause shifts in the radial positions of images 
near the ring radius. In particular, varying A^ and while leaving their sum unchanged affects 
the radial image displacements, but has little effect on the tangential image displacements (see 
Table |3|). The radial image shifts caused by these multipole terms depend on angular position, 
whereas those caused by do not. However, for lenses such as MG J0414+0534, in which there 
are only images at a few angular locations, the effects of and of A^ or B^ may compete, with 
large changes in /g compensating for small changes in A^ or B^. For systems with arcs or rings 
of lensed emission, information is available from a broader range of angles. The optical arc visible 
in MG J0414+0534 ( [Falco, Lehar fc Shapiro 19971 ) may further constrain MG J0414+0534's 



angular multipole moments, especially if its location can be measured with the same precision as 
the VLBA measurements used in this paper. 

That the radial profile parameter, /g, is so difficult to determine is unfortunate. It would 
provide information on how the angularly-averaged surface mass density decreases with radius 
near the Einstein ring radius, which could be used help choose a model value for (To(6). The 
quantity (To(&) is not directly constrainable from lensing, but it does affects the predicted time 
delays between images, a topic discussed in the next section. 



5.2. Implications for the time delays 

Each of the multiple images formed by a gravitational lens represents the source object at a 
different moment in its history. This is because the propagation time from source to observer is 
different for each image, due to the different path lengths and Shapiro delays experienced by the 
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ray bundles composing each image. In particular, the time delay by which an image at Fq, lags 
that at is, as computed by Narayan & Bartelmann (1996), 

Ai/3„ = t(r,) - t(r^) = (1 + ^l)^^ (^|Vr.<I>(r„)l2 - i| Vr,<I>(r^)P - <I>(r,) + $(r^)) (62) 

It is convenient to define a "dimensionless time delay" Ar^^ which depends only on the modeled 
lens potential and requires no values for redshifts or cosmological parameters, 

At0^ = l\VrMra)\^ - ^|Vr,$(r;3)P - <l>(r„) + cl>(r^) (63) 

and which is related to the time delay by a conversion factor depending on the redshifts and 
cosmology, 

At^o. = (1 + zi)^^^ATp^. (64) 

A measurement of the time delay between the flux variations in corresponding images and the 
lens redshift zi, when combined with a model that predicts the dimensionless time delay, thereby 
amounts to a measurement of the combination of angular diameter distances D^Ds / Dig. Since 
the relation between angular diameter distance and redshift depends on the values of i^o, f^mi and 
JIaj these cosmological parameters can be thereby constrained. The appeal of this well-known 



cosmological probe is that it does not rely on the usual intermediate distance indicators (Refsdal 
1964, 196q ). 



One problem with this idea arises from the mass-sheet degeneracy, which was discussed in 
section |3.4] . When the potential is transformed by the mass sheet degeneracy, the time delays 
between components are likewise transformed. Using the model potential of equation the 
model dimensionless time delay, Ar', is related to the true dimensionless time delay, Ar, by 

= 2 f 1 ^"""^ ^^^^ 

where o"o(6) is the angularly-averaged surface mass density at the Einstein ring radius. Unless the 
mass-sheet degeneracy can be resolved by determining o"o(6)/cTc in some independent fashion the 
time delays cannot be predicted, although the ratios of time delays between different image pairs 
can still be predicted. 

Although MG J0414+0534 has been extensively monitored at radio wavelengths, no flux 
variations have been observed that are large enough to permit an accurate time delay measurement 



( Moore fc Hewitt 1997 ). Nevertheless, this does not preclude radio or optical detections of time 
delays in the future, so it is important to understand how our models of MG J0414+0534 constrain 
the time delays. As discussed in the previous section, the best-fit model was 6 + A2 + B2 + A3 + A4. 
We used this model to make our single best prediction for the time delays of MG J0414+0534, by 
computing the dimensionless time delays between the 4 images of the brightest component (p). 
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The uncertainty in this particular model's predicted time delay due to the uncertainty 
in the measured image positions was estimated in the following manner. We re-computed 
the time delay with each parameter (one at a time) adjusted to its maximum and minimum 
values allowed by the stat- width confidence limits, while minimizing the over all the other 
parameters. The range between the highest and lowest time delays that were achieved during 
these parameter-by-parameter adjustments is our estimate for the uncertainty in the time delay 
for that particular model. 

However, we must also take into account the larger uncertainty in the predicted time delay 
caused by the uncertain choice of model. To estimate this uncertainty, we computed dimensionless 
time delays for a large subset of the models that were discussed in the previous section. For the 
11-parameter and 12-parameter models, we included models which adequately fit the observations 
when using the stat- width error estimates {N^^ < 3). The +A3 + M|"™ + and +A3 + B4 + 
models were excluded because, for these models, the entire confidence range for /g lies above 
/g > 2, well within the unphysical region /g > 1. Likewise, the +A3 + G22 model was not included 



because of its unphysically large value of IG22I, as discussed in section 4.2. The results are shown 
in Table |9|. Since these models make somewhat diff^erent predictions, any attempt to make a single 
prediction for the time delays must incorporate the uncertainty associated with the selection of a 
single model. Thus we have enlarged the error spread in our best predictions for the time delays 
to include the whole range of time delays predicted by all these models. The resulting predictions 
are: 
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(68) 



Here the "formal errors" represent the uncertainty in the parameters of the best-fit model, 
the "A1-A2 difference" is half the time delay between images Al and A2 (since the joint A1-A2 
light curve would probably be used to measure time delays), and "which model" refers to the 
uncertainty due to choice of model. The factor 2 (1 — <To(6)/crc) represents the uncertainty in the 
time delay due to the mass-sheet degeneracy, which can only be relieved by obtaining a reliable 
value of (TQ{b)/ac from other observational or theoretical sources. For a potential with the radial 
profile of an isothermal sphere, 2 (1 — cjo(^)/o'c) = 1- 
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To express the time delay as an actual number of days, the conversion factor of equation 
must be used. For MG J0414+0534, which has zi = 0.9584 and Zg = 2.639, this conversion factor 
takes the value At/Ar = 6.794 x 10^^ h^^ days, assuming the universe has 0^ = 1, = 0, and 
Hq = 75/175 km/s/Mpc. The values of this conversion factor for some other choices of cosmological 



parameters are tabulated in Table IC. In all cases, the "filled-beam" approximation was used 
to compute the conversion factor, in which the universe is assumed to have a perfectly smooth 
distribution of matter. The presence of dumpiness would require the angular-diameter distances 
to be re-computed (see e.g. Fukugita et al. 19921) . 



A promising way to reduce the "which model" uncertainty is to measure the time delay 



ratio AIab/ ^tAC- Predictions for this ratio are presented in Table 11 for various models. If this 
quantity is in accordance with the prediction of our best-fit model b + A2 -|- B2 -|- A3 -|- A4, and it 
can be measured to within 3%, then it would exclude all the other models listed and the "which 
model" error would fall away. Even if the ratio could only be measured to within 18%, it would 
exclude all but one other model, which would shrink the "which model" uncertainty in AIab to 
only 3.7% (from 62%). In this scenario, and using the conversion factor (eq. for = 1, 
JIa = 0, we would predict 
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6. Conclusions 

Upon first seeing the rich sub-structure in each VLBI image of MG J0414-I-0534, we were 
hopeful that such a large body of precise constraints on the modeling potential would lead to 
tight constraints on the mass distribution of the deflector and the predicted time delays between 
images. We hoped that the mathematical generality of the multipole- Taylor expansion (with slight 
modifications) would allow us to draw such conclusions without contamination from (perhaps 
faulty) astrophysical preconceptions. These hopes were only partly fulfilled. 

Once the best-fit parameters in the expansion are determined, it is difficult to know from 
which astrophysical source they arise. For example, our best-fit model b + A2 -|- B2 + A3 -|- A4 
seems to imply that either the mass distribution in the lensing galaxy is asymmetric and somewhat 
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quadrangular (boxy), or else that an unseen external perturber is partly responsible for the light 
deflection (as discussed in section It is possible, however, that the values of the best-fit 

model parameters may guide future interpretations of observations for this system, by indicating 
the possible directions of perturbing masses. 

Predictions of time delays, while quite well-constrained for any particular model potential, are 
limited by the uncertainty in selecting one of several viable models. In other words, the predictions 
are not limited by the precision of the positional constraints, but rather by the ability to satisfy 
those constraints with several alternative truncations of the multipole- Taylor series. In the case 
of MG J0414-I-0534, we found that one way to relieve this crucial source of systematic error is to 
measure ratios of time delays, which are predicted to have different values by different models. 

Despite these limitations, which afflict all lens modeling techniques to date, the multipole- 
Taylor expansion does indeed seem to be an appropriate form for multiple-image gravitational 
lenses such as MG J0414-I-0534. The simplest three-term truncation reproduces the observed 
image configuration far better than previously-used simpler models, and successively higher 
terms improve the fit by incrementally smaller amounts. The radial profile parameter could 
not be constrained, but this is likely to be the case for any point-modeling scheme applied to 
quadruple-image lenses. 

We believe the multipole- Taylor expansion could be usefully applied to other lenses in which 
the constraints occur at locations close to the Einstein radius, and the angular variation of 
the potential is expected to be fairly smooth. One serious problem with the efforts to date to 
determine Hq and other cosmological parameters by measuring time delays is that each known 
gravitational lens has usually been modeled in an individual and idiosyncratic manner. The 
multipole- Taylor expansion is one candidate for a very general modeling technique that could be 
applied to all the time-delay lenses, so that the results of these efforts could be sensibly combined. 

We thank Paul Schechter for many valuable discussions concerning lens modeling in general 
and this analysis of MG J0414-I-0534 in particular. J. N. Winn thanks the Fannie and John 
Hertz Foundation for financial support. This work was supported by grant AST96-17028 from the 
National Science Foundation. 
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Fig. 1. — High-resolution maps of the four images of MG J0414+0534, observed at 5 GHz with the 
VLBA, are shown superimposed on the 22 GHz VLA map of Katz, Moore fc Hewitt 1997| (kindly 
supplied by C. Katz). 
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Fig. 2. — High-resolution maps of the four images of MG J0414-I-0534, observed at 5 GHz with 
the VLBA, are shown with the four components p, q, r and s identified. The positions and fluxes 
for these components are reported in Table ||. 
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Fig. 3. — Observed and model positions for the images of MG J0414+0534. Clockwise from lower 
left, the panels show the vicinity of the images Al, A2, B and C. The model positions are plotted as 
small crosses (x). The error ellipses represent the 90%-confidence intervals for the observed image 
positions. The large error ellipses were computed using the stat-width error estimates, whereas the 
small error ellipses (most of which appear pointlike on a plot of this scale) were computed using 
the statistical error estimates. The different error estimates are described in |4.l| . 
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Fig. 4. — Model & + A2 + B2 + A3, fitted to MG J0414+0534. Tlie model parameters are illustrated 
graphically. North is up and east is left. An x marks the model deflector center of mass {gx,gy), and a 
ring is drawn at the model Einstein ring radius b. The error range on the ring radius is shown by additional 
rings drawn at the upper and lower confidence limits. (At this plot scale, neither the error rings nor the 
error range for the model galaxy position can be resolved.) Each m > 1 multipole moment is illustrated by 
arrows pointing in the possible directions to mass excesses that would cause such a multipole moment. The 
arrows are labeled with the number m of the multipole moment, and their lengths are proportional to the 
amplitudes of the multipole moments. Confidence ranges on the angles are indicated by small arcs at the 
tip of each arrow. Also shown are the R and /-band center positions of the defiector galaxy (+ marks), the 
position of Object X, the direction to a nearby group of galaxies (arrow), and the orientation of the galaxy 
isophotes (arcs near the ring radius), all from Falco, Lehar & Shapiro (1997) The locations of the observed 
and modeled components are also shown but cannot be distinguished at this plot scale. See Figure ^ for a 
closer view. The model sources are plotted northeast of the galaxy center. 
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Fig. 5— Model 6 + A2 + B2 + A3 + A4, fitted to MG J0414+0534. The model parameters are 
illustrated graphically, using the same conventions as in Figure |^. 
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Tabic 1: Centroid locations and flux densities of the elliptical Gaussian fits to the components of 
MG J0414±0534, derived from the VLBA observations described by Trotter (1998). The positions 
are relative to the correlation and fringe-fitting center at Al. The x-coordinate increases to the 
west and the y-coordinate increases to the north. The errors on the centroid positions are standard 
deviations on the x- and y-positions due to the thermal noise in the map, and the corresponding 
correlation between the x- and y-position errors. The flux errors are standard deviations due to 
the thermal noise in the map. These error bars do not include the effect of deconvolution error. 
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Einstein ring radius. Sensitive to total mass within ring. 


m 


= 


t = 2 




Sets surface mass density at ring, cro{b) (unconstraina- 
ble due to mass-sheet degeneracy). 
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Depends on ao{r) and its derivatives near r = b. For 
t = 3, gives falloff of surface mass density with radius, 
evaluated at the ring radius. See Eq. 44. 
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Dipole moment of mass exterior to ring (unconstraina- 
ble due to prismatic degeneracy). 
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exterior 




Sensitive to multipole moment of mass exterior to ring. 
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Sets COM of mass interior to ring. (The COM may be 
placed at the origin by setting Bi = 0.) 
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Sensitive to multipole moment of mass interior to ring. 
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Depends on (Jm{f) and its derivatives near r = b. 


m 


> 2 


t = 


TVTSUm 


M^™ = Am + Bm, sum of external and internal mass 



multipoles. (This parameter is from the parameteriza- 
tion in equation 7, rather than equation 47). 



Table 2: Physical significance of the parameters and terms in the modified multipole- Taylor series, 
equation 47. The first column lists the multipole moment m. The second column describes the 
radial dependence (the exponent t of p*, or whether the term is sensitive to mass interior or exterior 
to r = 5). The third column lists the name of the parameter. The last column describes the physical 
significance of the term. The parameters fo, /2, Ai, and Bi cannot be constrained from lensing, as 
explained in section 3.4. We set these parameters equal to zero during the model-fitting procedure. 
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+^ {-iflTT^} sinm0 - cosm^ 


Small, predominantly 
radial 


m 


> 1 


t = 


^ 1 ^ } {M^"<f' sin me - M^^' cos m^} 


Tangential, due to over- 








all strength of multipole 

m-"' = a;„ + b;„ 


m 


> 1 


t = l 


f- {-m} |M^f ' cos m6' + M^f sin 7716*} 


Predominantly 








+^ { } { Af^r sin me - M^f cos m^} 


radial. Affected by bal- 








ance between external 

and internal multipole 

lyrdiflf' _ A ' — R' 
m m m 



Table 3: The first derivatives of the terms in oiu parameterization of the lens potential. These are 
needed to calculate the displacement of the lens images from the source positions (eq. 3). The upper 
section of the table lists the five functional forms needed to represent the terms in the modified 
multipole- Taylor expansion of the potential (eq. 47) . The lower section of the table lists functional 
forms for terms in the original expansion (eq. 7). The identifications m and t are as in Table 2. 
The model parameters are written in Cartesian form, e.g. = xA^^, -|- yA'^. 
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m 


t 




m = 


t = 





m = 


t = i 


flfl i 1 I 


m = 


t > 2 


(t-i)! (1+p) /•!/*> 


m > 1 


exterior 


rr |-m(m - 1)(1 + | {^^c cos m0 + sin m6'} 
-r"" ^ TnyfTii — iji^i-|-pj^ -[/ij^j, cos TTiu -p -^yns sm Tnu J 
+{re + 0f) {m(m - 1)(1 + p)'-'^''^^} {A^,^ sin m6» - A^^^ cos m6l} 


m > 1 


interior 


rr { (i+^^)(i+2) } {B'mc cos m0 + B'^ sm m^} 
+00 [ } {S;^c cos me + B;,, sin me} 
+{re + er) {-^gg^} {B;,, sin me - S^^^cosm^} 


m > 1 


t>2 


{ ft-2)! } {^'mtc COS me + sm me} 
{ H - ("^' - *)^) } {G'mtc COS me + sin me} 
+(r0 + 0f) {- 1 + - 1)P)} {G;„tc sin me - G'^^, cos me} 


m > 1 


t = 


{ (if^ } {M--' cos me + M--' sin me} 
+(r0 + 0f){ j^fy] {M^^' smm9 M^^' cos m9} 


m > 1 


t=l 


00 [ (ify(l (m^ l)p)} |M^f cos me + M^f sin me} 
+{re + 0r) { (j^} {M;^f sin me - M^f cos me} 



Table 4: The second derivatives of the terms in our parameterization of the lens potential. These 
are needed for calculation of the magnification matrix. The notation is that of Table 3. The 
magnitude of the effect of each term on the inverse magnification matrix (eq. 55), may be seen by 
considering the exponent of p. 
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Model ^ l-K=^_lIII J[ 



2( 



point mass 


6^ Inr 





1/2 


1 


singular isothermal sphere 


br 


1/2 


1 





mass sheet with a = ac 


2' 


1 


oo 


-1 


power law 










a — > 0, point mass 


fe2 /r\2a 
2a U/ 








a = 1/2, isothermal 


a 


1 

2(l-a) 


l-2a 


a ^ 1, mass sheet 








power law with core radius 






1 s2+62 


(l-2a)6^-3s^ 


2a Vs^+''^/ 




2(l-a) 62 


s2+62 



Table 5: Comparison of several physically motivated monopole profiles. The first column is the 
lens potential. The second column is the surface mass density at the ring radius. If a lens with 
such a profile were fitted with the multipole- Taylor model with f!^ set to zero (which was the choice 
for /2 we employed) then the resulting scale factor 1 — k is listed in the third column. The fourth 
column lists the value of /g, the lowest-order monopole parameter that is sensitive to the radial 
mass profile. Note that if > 1, the surface mass density increases with radius near the ring 
radius. 
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Method for estimating positional errors 




Model 


N 


Ndof 


fit-size 
(an upper limit) 


statistical 
(a lower limit) 


Stat- width 
X' 


SIEP 

SIS+XS (6 + A2) 
PM+XS 


5 
5 
5 


19 
19 
19 


2.357e4 
1.680e4 
1.448e4 


5.402c3 
3.851e3 
3.317e3 


1.297e8 
1.078e8 
1.012e8 


2.976c7 
2.474e7 
2.321e7 


3.667o5 
1.910e5 
2.080e5 


8.412e4 
4.382e4 
4.771e4 




fe + A2 + B2 


i 


17 


2. 18502 


r).()i:)('i 


i.():i()(-r) 


:!.9()7el 


:!.!)!)()c:i 


!).():)lo2 





+B3 


9 


15 


1.356c2 


3.113cl 


1.143c5 


2.950c4 


2.022c3 


5.181c2 


+G22 


9 


15 


1.757c2 


4.149cl 


9.875c4 


2.549c4 


1.357c3 


3.464c2 


+/3 


8 


16 


1.108c2 


2.371cl 


1.040c5 


2.601c4 


1.438c3 


3.554c2 


+G12 


9 


15 


1.400c2 


3.227el 


1.430c5 


3.691c4 


8.457c2 


2.145e2 


+A4 


9 


15 


1.011c2 


2.224cl 


7.958c4 


2.054c4 


7.493c2 


1.896c2 




9 


15 


1.235c2 


2.801cl 


8.522c4 


2.200c4 


3.913c2 


9.717cl 


+B4 


9 


15 


4.168cl 


6.889 


5.551c4 


1.433c4 


1.983c2 


4.733cl 


+A3 


9 


15 


1.619cl 


0.308 


1.831c4 


4.723c3 


1.351c2 


S.lOOcl 




9 


15 


8.783 


-1.605 


1.603c4 


4.135c3 


7.038cl 


1.430el 



-I-A4 + B4 


11 


13 


1.917el 


1.712 


3.370e4 


9.344e3 


1.297e2 


3.235el 


-l-MI™ + /a 


10 


14 


8.113 


-1.573 


1.596e4 


4.261e3 


6.284el 


1.305el 


-t-MI""" + B4 


11 


13 


5.087 


-2.195 


6.435e3 


1.781e3 


5.647el 


1.206el 


+M|""" + MT'" 


11 


13 


3.795 


-2.553 


6.130e3 


1.697e3 


2.862el 


4.331 


-hMI"'" + A4 


11 


13 


2.459 


-2.924 


3.998e3 


1.105e3 


1.239el 


-0.1692 


+Mf - + G22 


11 


13 


3.285 


-2.694 


8.421e3 


2.332e3 


1.179el 


-0.3348 


+M|"'" + G12 


11 


13 


2.527 


-2.905 


7.794e3 


2.158c3 


1.145cl 


-0.4291 


+A3 + Ml™ 


11 


13 


9.435 


-0.989 


1.223c4 


3.389c3 


5.506cl 


1.167cl 


+A3 + B4 


11 


13 


7.374 


-1.560 


1.050c4 


2.908c3 


4.732el 


9.519 


+A3 + /3 


10 


14 


9.144 


-1.298 


1.631c4 


4.356c3 


4.183cl 


7.439 


+A3 + G12 


11 


13 


7.254 


-1.594 


1.528e4 


4.235c3 


3.612cl 


6.413 


+A3 + Ba 


11 


13 


1.654 


-3.147 


4.176e3 


1.155c3 


8.960 


-1.120 


+A3 + G22 


11 


13 


4.783 


-2.279 


5.986c3 


1.657c3 


8.387 


-1.279 


+A3 + A4 


11 


13 


1.257 


-3.257 


3.863c3 


1.068c3 


6.073 


-1.921 




-I-A4 + B4 + /a 


12 


12 


1.814el 


1.773 


2.871e4 


8.285e3 


8.107el 


1.994el 


-l-MI™ + B4 + /a 


12 


12 


4.889 


-2.053 


6.342e3 


1.827e3 


3.179el 


5.711 


+M|""" + Ml""^ + fs 


12 


12 


3.791 


-2.370 


5.948e3 


1.714e3 


2.455el 


3.622 


+Mr'^ + A4 + /3 


12 


12 


2.321 


-2.794 


3.798e3 


1.093e3 


1.234el 


0.099 


+Mr'^ + G22 + /a 


12 


12 


1.364 


-3.070 


2.795e3 


8.033e2 


9.373 


-0.758 


-FMf" + G12 + fs 


12 


12 


2.203 


-2.828 


5.505e3 


1.586e3 


8.495 


-1.012 



-^Aa 


+ G12 + /a 


12 


12 


6.160 


-1.686 


1.445c4 


4.169c3 


3.431cl 


6.439 


+A3 


+ Ml™ + /a 


12 


12 


3.376 


-2.490 


7.549c3 


2.176c3 


1.696cl 


1.431 


+A3 


+ Ba /a 


12 


12 


1.571 


-3.011 


3.989c3 


1.148c3 


8.193 


-1.099 


+A3 


+ G22 + /a 


12 


12 


3.163 


-2.551 


4.895c3 


1.410c3 


7.334 


-1.347 


-HA3 


+ B4 + /a 


12 


12 


2.000 


-2.887 


4.914c3 


1.415c3 


6.157 


-1.687 


H-A3 


+ A4 -f /3 


12 


12 


0.720 


-3.256 


2.306c3 


6.622c2 


2.721 


-2.679 



Tabic 6: Results of fitting multipole- Taylor models to MG J0414-I-0534. For each model, the 
number of parameters N and number of degrees of freedom N^qf are hsted, along with the 
minimum value of obtained and the number of standard deviations N^j separating the model 
and constraints. Models with negative N^r oversatisfy the positional constraints. Three different 
estimates for the positional errors were used and are reported separately (as discussed in section 4.1). 
The models are named as described in section 4.2. 
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denector positions: 








(W and N of the correlation 


9x = 


1.0788 ± 0.0020 


arcseconds 


center at Al) 


9y = 


0.6635 ±0.0012 


arcseconds 


ring radius: 


ft - 1 1 474+U.UU25 

— J---"-^' ^-0.0026 


arcseconds 


internal quadrupole: 


B'2 = 


U.Ui04Z_Q QQQgg 




B2 = S2(x cos 2-052 +ysin2'0B2) 




— U. / lo_Q Q28 


radians N of W 






= 69.57tO|J 


degrees E of N 


external quadrupole: 


A' — 


0.04478^|j:S 




A2 = A2 (x cos 2ipA2 + y sin 2ipA2 ) 




= -o.5i3ini^ 


radians N of W 






= 75.29^0:1 


degrees E of N 


external m = 3 multipole: 




o.oio8ot|]:S 




A'3 = A'^ (x cos 3-0^3 + y sin Sip A3 ) 




= 2.678;°;°i 


radians N of W 








degrees E of N 



external m = 4 multipole: A'^ = 0.00415+^;^^^^^ 

= >l^(x cos 4^A4 + y sin 4^a4 ) ^i'A^ = 2.020lg:|]^^ radians N of W 

V'A4 = 28.93l;J;^^ degrees E of N 
45° + ipAi= 73.931^;^^ degrees E of N 



undeflected source positions: 


Sp 


X = 


0.8795 


arcseconds 




Sp 


y = 


0.7549 


arcseconds 




Sg 


X = 


0.8815 


arcseconds 




Sg 


y = 


0.7559 


arcseconds 






X = 


0.8928 


arcseconds 






y = 


0.7626 


arcseconds 




Ss 


X = 


0.8687 


arcseconds 






y = 


0.7509 


arcseconds 



Table 7: Best-fit model parameters for the deflector model +A3 + A4. See Figure 5 for a 
graphical depiction of these model parameters. The limits provided for each parameter are formal 
68.3% confidence limits, computed as described in section 4.1. For completeness, the lower half 
of the table lists the undeflected source positions for the model sources; these are a by-product 
of modeling the deflector gravitational potential, and no confidence ranges on the source positions 
were calculated. The image magnifications predicted by this best-fit model are listed in Table 8. 
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Image Component Magnification 



Al 


P 


12.4 


Al 


q 


12.6 


Al 


r 


15.3 


Al 


s 


12.2 


A2 


P 


-14.1 


A2 


q 


-14.2 


A2 


r 


-16.4 


A2 


s 


-14.5 


B 


P 


5.14 


B 


q 


5.06 


B 


r 


4.66 


B 


s 


5.58 


C 


P 


-1.84 


C 


q 


-1.85 


C 


r 


-1.89 


C 


s 


-1.78 



Table 8: Image magnifications, as predicted by the best-fit model -I-A3 + A4, which is specified 
by the parameters in Table 7. 



-43- 



Model 


^'^'a\vA2v ^ ^^^^ 




ArAi,cp X 10'' 




+A3 


u.ut)iy_Q 0009 
n nc^88+o ooo9 

U.Ul3»»_o.0009 


9 Q9K+0.038 
^.ozu_o 039 

/loi +0.037 


8.661°:°^ 
9.96l°:li 


10.991°:!^ 

12.44l°:l« 


+M|"'^ + A4 
+Mf + G22 
+Mf + G12 


0.0660t°:°°l^ 
0.06191°:°°!^ 
0.0485lHgii 


z.oiu_o 047 

Z.DiU_o.077 
-2 061+° °59 


9.87l°:i 

11 oQ+0.20 
il.JJ_0 17 

8.38l°:i 


12.49l°:i 
13.94l°:le^ 

io.44l°:i 


+A3 + B3 
+A3 + G22 


0.05891°:°°!? 
n^sfi+°-°°°3 

U.UOOD_o.0004 


-2.4»b_o.028 
l.yoO_o 020 


10.691°:^° 

1 1+0.15 
^•-^■'■-0.15 


13.181°:- 
io.iol°:l? 


+A3 + A4 


n 051 1^+0 0019 


1 ono+0.034 
i.OUZ_o 036 


10.44+gJ^ 


12 25+°-l'^ 


+Mf - + A4 + /a 


0.0680ini?i 


-2.694l°:?ol 


io.iol§:i 


12.80l?-J^ 


+Mf - + G22 + /3 
+Mf- + G12 + /3 


0.0709l0-0^3i 

0525+°-°°^^ 
u.uozo_o 0020 


-2 qi 9+0-664 
z.yiz_2.092 
9 909+0-075 
Z.ZUZ_Q 044 


8.64l^:^« 
8.981°:^? 


11 ^f^+7'57 
n 10+0.21 

-•--^•-^0-0.38 


+A3 + B3 + /3 

+A3 + G22 + /a 


0.0493l°:°J^} 
0.05881°:°°^! 


9 -|Oq+0.161 

Z.100_Q 4]^0 
9 1 QQ+0.194 
Z.iOO_o.l74 


9.49l-« 

8 ^7+0-63 
-0.63 


il.DZ_j^ 08 

10.75l°:?9^ 


+A3 + A4 + /3 


0.04701°:°°?? 


1 70/I +0.029 
1. /O^_0,029 


q 1 Q+0-61 


10.86l°:i 



Table 9: Dimensionless time delays for models fitted to MG J0414+0534, using the stat-width 

estimates of the observational uncertainties. The dimensionless time delay is defined in equation 64. 
The results are reported for the brightest component, p. Variations would first be observed in image 
B, then Al, A2, and finally C. 



^rn M/^T X 10~ 

(^75' days) 
1 6.794 
0.3 0.7 7.654 
1 7.166 
8.785 



Table 10: Values of the conversion factor between the time delay and the dimensionless time delay, 
in various cosmological scenarios. 
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Model 


Ar BpAip/ Ar AipCp 


Ar BpA2p/ ^TA2pCp 


+A3 


0.268t°:™3 


n 27«+o.ooa 
u.z(u_o.ooa 




0.249^°:°°^ 


>J-ZOD_o 002 


+Mf + A4 


U-^OO-O.Oll 


n 974+0-011 
u-^' ^-0.011 


+Mf + G22 


n 2sn+°-°°'^ 

U.^OU_o,oo8 


r, 907+0.007 

U-^-J' -0.008 


+Mf + G12 


n O4fi+o-oo2 

U.Z41:D_o.0O2 


n 9t:q+0.002 
U.zoo_o,oo2 


+A3 + B3 


n 000+0.011 

^•^^"-•-0.011 


u.zoy_o 012 


+A3 + G22 


r, 24P;+0.004 
U.Z40_Q 004 


n 9c;q+0.004 
U.ZOO_o.004 


+A3 + A4 


n -, 70+0.006 


n 1 78+0.006 
U.i/S_0.006 


+Mf- + A4 + /3 


U-^O' -0.012 


n 97c;+0.015 
U.z/o_o 013 


+Mf- + G22 + /a 


PI 007+0.020 
u.oor_o.oi7 


0.348;^?^ 


+Mf- + G12 + /a 


U.Z40_o 001 


n ocq+O.OOl 
u.zdo_o,ooi 


+A3 + B3 + /3 


r, 9r,c-+o.oia 
U.^zo_o 008 


Q 201+0-014 
U.zoi_o,oo8 


+A3 + G22 + /a 


Q 2c;c;+0.016 
U.ZOO_o 014 


U.ZOO_o 015 


+A3 + A4 + /a 


u.iyu_o 012 


n 196+0012 

U.iyD_o 012 



Table 11: Time delay ratios for models fitted to MG J0414+0534. The quoted errors were 
computed using the stat- width estimates for the positional errors. 



